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ABSTRACT 

In calculating the position vector of the Moon in on-board flight software, 
one often begins by using a series expansion to calculate the ecliptic latitude 
and longitude of the Moon, referred to the mean ecliptic and equinox of date. 
One then performs a reduction for precession, followed by a rotation of the 
position vector from the ecliptic plane to the equator, and a transformation 
from spherical to Cartesian coordinates before Anally arriving at the desired 
result: equatorial J2000 Cartesian components of the lunar position vector. An 
alternative method is developed here in which the equatorial J2000 Cartesian 
components of the lunar position vector are calculated directly by a series ex- 
pansion, saving valuable on-board computer resources. 

INTRODUCTION 

The calculation of the orbit of the Moon is one of the oldest problems in 
celestial mechanics. Its solution has had great historical significance as a test 
of Newton’s theory of gravity, with much of the early work on the problem 
having been done by Newton himself in his discussion of the two- and three- 
body problems in Book I of the Princ.ipia. In past centuries, accurate predictions 
of the position of the Moon have also been of great practical interest as a 
navigational aid for seafaring vessels, prompting the English government and 
scientific societies to offer rewards for accurate lunar prediction tables. 1 The 
resulting body of work developed during the eighteenth and nineteenth centuries 
forms the basis of the lunar theory still in use today. 

Modern lunar theory was first developed bv G.W. Hill 2-0 in 1878, and later 
expanded and improved by E.W. Brown 6 in 1890. The problem of lunar motion 
addressed by Hill and Brown is a surprisingly difficult one; while the underlying 
physical laws are very simple, the motion itself is quite complex.' -11 The basic 
motion of the Moon around Earth is affected by many strong perturbations such 
as those due to the Sun, the other planets, and Earth's equatorial bulge. These 
perturbations result in an advancement of the line of apsides of the lunar orbit, 
a regression of the line of nodes, and other periodic perturbations superimposed 
on these motions. For high accuracy, it is necessary to compute' hundreds of 
periodic variations in the motion, although computing only the most important 
few terms results in a level of accuracy that is adequate for flight, software use. 

Then' have been two major reasons for calculating the position of the Moon 
in spacecraft oti-board computer flight software. First, one often wishes to write 



flight software to prevent the spacecraft from pointing sensitive instruments at 
the Moon, which can have an apparent magnitude as bright as —12 at full 
Moon. 12 Second, one may require the flight software to calculate stellar aberra- 
tion corrections. 13 For high accuracy, this requires calculating the velocity vector 
of Earth with respect to the Earth-Moon barycenter. which in turn requires a 
calculation of the lunar velocity vector. If the flight software can calculate a 
lunar position vector, then this velocity vector may be found bv differentiating 
the lunar position vector with respect to time. 

REVIEW OF CURRENT MODELS 

A number of approaches for calculating a lunar position vector are currently 
used by spacecraft flight software. In the flight software for the Hubble Space 
Telescope’s DF-224 flight computer, for example, one finds the position of the 
Moon using a simple two-bodv model. The standard two-bodv calculations 14 
are modified somewhat to allow for the motion of the nodes and apsides of the 
lunar orbit. A new set of orbital elements is uplinked from the ground every 
few days to keep the error in the model to within acceptable limits, on the order 
of 1°. While this model is not highly accurate, it has the virtue of being very 
fast -a necessity for the 1970s-vintagc flight computer. 

An approach commonly used with more modern flight computers is based on 
the low-precision formulae given in the Astronomical Almanac. 10 ' 16 This model 
is based on earlier w'ork done by the Almanac Offices of the United States and 
United Kingdom 17 and by Eckert, Walker, and Eckert, 18 all of which are based 
on Brown’s lunar theory. 6 In this model, one begins by using series expansions 
to calculate the ecliptic longitude A, ecliptic latitude' T and horizontal parallax 
7 r of the Moon, referred to the mean ecliptic and equinox of date: 

A = 218?32 4 481 2G7?883 t 

4G?29sin(477 198?85 t 4 134%) 

— l?27sin(— 413 335?38 t 4 259°2) 

4() o GGsin(890 534°23 t 4 235?7) 

+0?21 sin(954 397?70 /. 4 - 2G9°9) 

— {J?19sin(35 999°05 t 4 357?5) 

— ()?1 1 sin(9GG 404?()r> t 4 18G°G) , (1) 

fi = +5?13 sin (483 202°()3 t 4 93°3) 

4()?28sin(9G0 400?87 t 4 228?2) 

— 0?28sin(G 003°18 t 4 318°3) 

— ()°1 7 sin( -407 332?20 t 4 217%) , 
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7 r = ()?0. r j()8 

+0?()518cos(477 198?85 t + 134°9) 

+0°<m'3<:os(-413 335?38 t + 2, r )!)?2) 

+n?0078 (x>s(890 534 °23 t + 233?7) 

+0?0028 cosfOM 397?70 t + 209°!)) . (3) 

Tho horizontal parallax 7r {■ivcs the Earth-Moon distance r: 

r = ■ (4) 

Sill 7T 

where 7?® = 0378.140 km is the equatorial radius of Earth (IAU 1976 value). 19 

Having found the lunar ecliptic mean-of-date coordinates, one must then 
perform a reduction for precession to epoch J2000 (2000 January 01 12:00:00 
Barycentric Dynamical Time) to find the ecliptic J2000 coordinates (Ao, /7 o). 


To sufficient precision, this may he found using the formulae -0 

fj 0 = {j — hi sin (A -he), (5) 

A 0 = A - a 4- />eos( A 4- e) tan 1% , (6) 

where the precession constants a, b, and r are given by 

a - 1?396971 t -h ()?{)()() 3086 t 2 , (7) 

h - 0?013056 t - ()?()0() 0092 t 2 , (8) 

e = 5?123 62 - 1?155 358 t - 0°000 1964 t 2 , (9) 

and where t is the time in Julian centuries (cy) of 36 525 days from J2000: 

t = (JDE - 245 1545.0)/36 525 , (10) 


and JDE is the ephemeris Julian day. 

The remaining step is to rotate the coordinates from the plane of the mean 
ecliptic of J2000 to the mean equator of J2000. and to convert from spherical 
polar to Cartesian coordinates: 


r cos /^o cos A o . 

(ii) 

r (cos /Jo sin Aq cos 5 q — si 11 fto s i n 5 o) ? 

(12) 

r (cos iio sin Ao sin zq 4- sin tio cos sq ) , 

(13) 


when' r is given hv Eq. (4) and 6 0 = 23° 26' 21'.'448 is the obliquity of the 
ecliptic at J2000 (IAU 1976 value). 21 

This model has very good precision for on-board flight software use: the mis 
error in the lunar position is about 0?1 1, with a maximum error of about 0?35. 



A NEW MODEL 


Many of the equations involved in computing the position of the Moon us- 
ing - the method just described involve what is essentially a coordinate transfor- 
mation, from ecliptic mean-of-date coordinates to equatorial ,12000 Cartesian 
coordinates. In this paper, I investigate the possibility of calculating the equa- 
torial J2000 Cartesian coordinates directly bv series expansions similar to Eqs. 
(1 3), thus eliminating the need for performing the coordinate transformations 
in on-board flight software. 

We begin by assuming that each of the J2000 equatorial Cartesian coordi- 
nates X n may be represented by Fourier sine series: 

JV„ 

X n ” ^ ^ a nm sm(ij nrn t H~ ^nm ) ? ( -Id) 

m — 1 

where X x = X, Xo = Y, and X 3 = Z; N n is the order of the series for X n . 
We now need to find the amplitudes a nm , frequencies w nm , and phase constants 
fi nrn . This may be done by fitting these parameters to the DE200 ephemoris 
model 22,23 using an exhaustive search. DE200 is an ephemeris model developed 
at the Jet Propulsion Laboratory, and has been used to produce tables in the 
Astronomical Almanac since 1984. It calculates Cartesian coordinates of Solar 
System objects, referred directly to the mean equator and equinox of J2000. 

For each coordinate, the terms of the series in Eq. (14) may be found one 
at a time by simultaneously fitting the parameters a nmi cj nm , and ^nm over a 
grid of possible values to the DE200 model. An algorithm for accomplishing 
this involves calculating the error c auJ s between the DE200 model and a "test 
model” asin (uft 4- 6) using each combination of parameters a , and 

for U — a.min to Clmax 

foi" (jJ — tO t ^rnax 

for = firnin tO fi-max 

faujS = S ^ 2000 [^^^ 200 (0 “ as\n(u;t -b ^5)]“ • 

where the summation is over 2 16 points covering the interval a.d 2000 2100. 
The smallest, error e auJ 6 found gives the best fit parameters a, u, and 6. This 
process may be repeated several times over successively smaller search ranges 
and finer grid spacings in order to find more significant digits for the parameters. 
Once a term has been found, it is subtracted from the DE200 data, and the whole 
process repeated on the remaining data to find the next term in the series. 

In the model given by Eq. (14), we assume that the amplitudes a nrn are all 
positive, so that amplitudes may bo searched over a grid of values between 0 and 
the maximum in the data set. The amplitudes may be assumed t,o be positive 
without loss of generality by allowing the phase constants fi nm to be searched 
over the entire range 0 to 27r: since - sin 0 = sin(0 4- tt), any potential minus 
sign in the? amplitude is simply absorbed as an extra tt radians added to the 
phase constant. 



Fourier Transform of Lunar X Coordinate 



Figure 1. Fourier spectrum of lunar X coordinate (a.d. 2000 2100). 


Determining a search range for the frequencies cj nm is somewhat more com- 
plicated than it is for the amplitudes and phase constants. A search range for 
cj nm may be determined by examining the peaks in the Fourier transform X n (u ;) 
of the DE200 data: tQO 

X n {u) = I X n {t)e tujt dt , (15) 

./ — oo 

where X n (t.) is the position coordinate at time and lj is the angular frequency. 
This Fourier transform may be calculated by using the DE200 model to compute 
the lunar position vector at N discrete time points t k , then finding the discrete 
Fourier transform X n (tJp): 

:\ r — 1 

Xn{* P ) = E X n (t. k ) . 

k = i) 


(io) 



where X n (tk) is the position vector at time point f*., cu p = 2 irp/ty is the angular 
frequency, and p = 0, 1, 2, . . . , N - 1. For this study, N — 2 14 time points were 
chosen over the time interval a.d. 2000 2100; the magnitude of the resulting 
Fourier transform |Xi(u; ; ,)| for X is shown in Figure 1. For each term in the 
series expansion (Eq. 14), a search range is taken around one of the peaks in 
the Fourier spectrum. 

This exhaustive' search process, which is essentially a curve fit to the DE200 
model, required about one week of computer time to find each term in a series, 
and some five months of computer time to find the complete solution to seven 
terms per series. The final results are: 

X = 383.0 sin (8399.685 t + 5.381) 

+ 31.5 sin (70.990 t + 6.169) 

+ 10.6 sin (16 728.377 t + 1.453) 

+ 6.2 sin (1185.622 t + 0.481) 

+ 3.2 sin (7143.070 t + 5.017) 

+ 2.3 sin (15613.745 t + 0.857) 

+ 0.8 sin (8467.203 t + 1.010) x 10 6 m , (17) 

Y = 351.0 sin (8399.087 t + 3.811) 

+ 28.9 sin (70.997 t + 4.590) 

-I- 13.7 sin (8433.400 t + 4.700) 

+ 9.7 sin (10 728.380 t + 0.105) 

+ 5.7 sin (1185.007 t + 5.104) 

+ 2.9 sin (7143.058 t + 0.300) 

+ 2.1 sin (15013.755 t. + 5.505) x 10 6 in . (18) 

Z = 153.2 sin (8399.072 t, + 3.807) 

+ 31.5 sin (8433.404 t. + 1.029) 

+ 12.5 sin (70.990 t + 4.595) 

+ 4.2 sin (10 728.304 t + 0.102) 

+ 2.5 sin (1185.045 t + 5.107) 

+ 3.0 sin (104.881 t + 2.555) 

+ 1.8 sin (8399.110 t + 0.248) x 10 6 in . (19) 

where all angles are given in radian s for eonvenienee of use in software, t is the 
time in Julian centuries from J2000 given by Eq. (10), and A. V. and Z are 
the Cartesian components of the lunar position vector, referred to the mean 
equator and equinox of J2000. The terms are arranged in order of decreasing 
contribution to the reduction in the error of the model. 



One of the primary advantages of this model is that it allows a lunar ophemeris 
to be programmed in flight software using very little rode. Using Eqs. (17 19). 
an entire lunar ophemeris model may be programmed in just a few lines of C 
code: 

for (n=0; n<3; n++) 

{ 

x [n] « 0.0; 

for (m=0; m<7 ; m++) 

x [n] += a [n] [m] *sin(w [n] [m] *t+delta [n] [m] ) ; 

> 

Calculations for the reduction for precession, rotation from the ecliptic to the 
equator, and transformation from spherical polar to Cartesian coordinates have 
essentially been “absorbed” into the series coefficients, and so do not need to 
be performed explicitly. 

DISCUSSION OF THE NEW MODEL 

An examination of the frequencies in the terms of the Astronomical Almanac 
model of Eqs. (1 3) and of the new model of Eqs. (17 19) gives some interesting 
insights into the lunar motion. The frequencies in the Astronomical Almanac 
model are all computed as functions of the mean anomalies and mean longitudes 
of the Sun and Moon, 16 while the frequencies in the model given by Eqs. (17 19) 
are determined entirely by a curve fit. We examine the origins of some of the 
more prominent frequencies in both models below. 

Anomalistic Month 

The dominant term in the expressions for the ecliptic longitude A (Eq. 1) 
and horizontal parallax tt (Eq. 3) have a frequency of 477198.85 deg e:y~U 
In deriving the Astronomical Almanac series, this frequency was computed as 
the rate of change of the Moon’s mean anomaly. Since the mean anomaly is 
measured in the plane of the orbit from the perigee point, one complete cycle 
of the mean anomaly requires the same amount of time as the Moon’s motion 
from its perigee point to its next perigee'. It comes as no surprise, then, that 
this frequency of 477198.85 deg cy -1 is equal to one revolution per anomalistic 
month of 27.554 550 days, where an anomalistic month is the time required for 
the Moon to move from perigee to perigee. 

Draconic Month 

For the ecliptic: latitude' (i (Eq. 2), the dominant term has a frequency of 
483 202.03 deg cv _l . This was computed as the rate of change' of the' Moon’s 
mean longitude', which is measuml from the vernal eepiinox to the ascending 
node along the ecliptic plane, then from the node' to the Moon along the' orbit 
plane'. The' Moon will have ji = 0 only wlie'ii it, is at one of the' nodes of the 
orbit, and it will next; have' ii = 0 again (crossing the' node in tlu' same' direction) 



when it returns to the same node again. We might therefore expert that the 
dominant term in the expression for the ecliptic latitude will he the time required 
for the Moon to move from an orbital node hack to the same node. Indeed, the 
frequency of 483 202.03 deg cv _l is equal to one revolution per draconic month 
of 27.212 221 days, where a draconic month is the time required for the Moon 
to move from an orbital node hack to the same noth'. 

Sidereal Month 

In the series for X , Y, and Z in the new model (Eqs. 17 19), on the other 
hand, the dominant terms all have a frequency of about 8399.085 rad cy -1 , 
which is equal to 1 revolution per sidereal rnonth of 27.321002 days, where a 
sidereal month is measured with respect to the fixed stars. This is a reflection of 
the model having its coordinate system fixed in space (mean of J2000 equatorial 
coordinates). 

Motion of the Apsides 

A comparison of the model of Eqs. (1 3) with the new model of Eqs. (17 
19) shows that the new model includes an important term that does not appear 
in the conventional model, having a frequency of about 70.99 rad cy -1 . This 
frequency reflects the motion of the line of apsides of the lunar orbit. The 
expected frequency of this motion may be computed from the periods of the 
anomalistic and sidereal months: 

2tt 2tt 

sidereal mo. anomalistic mo. 

( 27T 2n \ days 

_ x 30 525 

\ 27.321 M2 d 27.554 550 d ) cy 

= 70.9932 rad cy -1 (20) 

in close agreement with the frequencies found using the curve fit. 

ERROR ANALYSIS 

The results shown in Eqs. (17 19) have been checked against the DE200 
ephemeris model by using DE200 to generate lunar X , Y, and Z coordinates at 
2 20 (over one million) time points between a d 2000 January 1 and a.d. 2100 
January 1, corresponding to roughly one point every fifty minutes for 100 years. 
The model shown in Eqs. (17 19) was run at the same time points, and the 
results compared with the DE200 results. This error analysis shows an rrns 
position error between DE200 and the new mode*! of Eqs. (17 19) of 0?341, and 
a maximum error of 1?033. 



CONCLUSIONS 


Three lunar epliomeris models for on-board flight software use have been 
discussed. A modified two-bodv model is very fast, but is of low precision 
and requires constant maintenance in the form of periodic updates of orbital 
(dements from the ground. The model currently in common use. which is based 
on the low-precision formulae in the Astranouri('alAlrnana(\ is of very good 
precision and will run indefinitely without ground intervention, but requires code 
to convert the calculated ecliptic mean-of-date coordinate's to equatorial J2000 
Cartesian coordinates. The method developed in this paper is of intermediate 
precision, requires the least code of the three, and will also run indefinitely 
without ground intervention. It may have applications for small missions where 
computer resources are limited and its precision is acceptable. 
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